%% Preliminary Inputs
% Function Handle 
funcHandleOne=@modBMJ_newBaseDepU; 
funcHandleTwo=@modBMJ_newBaseDepU; 

addsol=[]; 
cucd = cd;

% Name and location of the Parameter File
mode.path=fullfile(cucd); 
mode.filename  = 'callParams';
mode.filename1 = 'callParams_01mu'; 
mode.filename2 = 'callParams_99mu'; 

% Number of horizons for the IRFS 
dim.horIRF=16; 

% Number of horizons for the forecast error decompositions 
dim.horDec=20; 

% Other Possibliities
DoIRF  = true;
DoCSum = false; 
closeplots = true;

% Customize the IRF Output
ShockstoPlot    = [1:14];
VartoPlot       = [1:27 41 44 53];
FG.numSignal    = 4;
PlotSignals     = false;


parPerturb.name='taupi'; 
parPerturb.vals=[1.25 3]; 


% Variables to cumsum
VarCumSum = [
            1  % Consumption
            2  % Hours
            5  % Investment
            18 % Wage
            %38 % GDP
            ];

% Names of Variables



home
disp(' ')
disp('------------------------------')
disp(' ')
%==========================================================================
%% Preliminaries
if closeplots
    close all
end

%% Load XLS File with Mode 
cd(mode.path);
eval(mode.filename); 
cd(cucd); 

%% Run Model
% Model Names
namesModels={[parPerturb.name,'=',sprintf('%2.4f',parPerturb.vals(1))]...
             [parPerturb.name,'=',sprintf('%2.4f',parPerturb.vals(2))]}; 
fldnm=parPerturb.name; 
paramBase.(fldnm)=parPerturb.vals(1);

[GGOne,RROne,CONSOne,euOne,SDXOne,ZZOne,ACont1,ALag1,ALead1,names]=feval(funcHandleOne,paramBase); 
if isequal(euOne,[1;1])==false; 
    error('Intedeterminacy'); 
end 

paramBase.(fldnm)=parPerturb.vals(2);
[GGTwo,RRTwo,CONSTwo,euTwo,SDXTwo,ZZTwo,ACont2,ALag2,ALead2]=feval(funcHandleTwo,paramBase); 
if isequal(euTwo,[1;1])==false; 
    error('Intedeterminacy'); 
end 


%% IRFs
[stdObsOne,stdStatesOne,irfObsOne,irfStatesOne]=...
    irfandthmom(GGOne,RROne,SDXOne,ZZOne,dim.horIRF,dim.horDec,zeros(size(ZZOne,1),1),0); 

[stdObsTwo,stdStatesTwo,irfObsTwo,irfStatesTwo]=...
    irfandthmom(GGTwo,RRTwo,SDXTwo,ZZTwo,dim.horIRF,dim.horDec,zeros(size(ZZTwo,1),1),0); 

if DoIRF
disp('Plotting IRFs...')

if any(~ismember(VartoPlot,1:size(irfStatesOne,1)))
    error(['Variables don''t exist in irfStatesOne, size = ' num2str(size(irfStatesOne,1))])
elseif any(~ismember(VartoPlot,1:size(irfStatesTwo,1)))
    error(['Variables don''t exist in irfStatesTwo, size = ' num2str(size(irfStatesTwo,1))])
elseif any(~ismember(ShockstoPlot,1:size(irfStatesOne,3)))
    error(['Shocks don''t exist in irfStatesOne, size = ' num2str(size(irfStatesOne,3))])
elseif any(~ismember(ShockstoPlot,1:size(irfStatesTwo,3)))
    error(['Shocks don''t exist in irfStatesTwo, size = ' num2str(size(irfStatesTwo,3))])
elseif any(~ismember(ShockstoPlot,1:size(namesPages,1)))
    error(['Shocks don''t exist in namesPages, size = ' num2str(size(namesPages,1))])
elseif any(~ismember(VartoPlot,1:size(namesRows,1)))
    error(['Variables don''t exist in namesRows, size = ' num2str(size(namesRows,1))])
end
fonts.axes=8; 
fonts.titfont=9; 
ghvec=plot_TwoIrfsTodd(irfStatesOne(VartoPlot,:,ShockstoPlot),irfStatesTwo(VartoPlot,:,ShockstoPlot),...
    names.states(VartoPlot),namesPages(ShockstoPlot),namesModels,'',fonts); 

disp(' ')
disp('------------------------------')
disp(' ')

ghvec=ghvec( (ghvec~=0) ); 
delete IrfEps*
for ii=1:length(ghvec);
    print(ghvec(ii),'-dpsc','IrfEps','-append');
end
%%
disp('Done!')
disp(' ')
disp('------------------------------')
disp(' ')


end

%% CumSum Plots
if DoCSum
    disp('Plotting Cumulative Graphs...')
    VarCumSum = VarCumSum';
    LevelOut = cell(2,1);
    for i = 1:2
        LevelOut{i} = zeros(numel(VarCumSum),dim.horIRF+1,numel(ShockstoPlot));
    end
    
    tempexp1 = exp(irfStatesOne);
    tempexp2 = exp(irfStatesTwo);
    for i = 1:numel(VarCumSum)
        for k = numel(ShockstoPlot)
            LevelOut{1}(i,:,k) = cumsum(tempexp1(VarCumSum(i),:,ShockstoPlot(k)));
            LevelOut{2}(i,:,k) = cumsum(tempexp2(VarCumSum(i),:,ShockstoPlot(k)));
        end
    end
    
    fonts.axes=9;
    fonts.titfont=10;
    ghvec=plot_TwoIrfs(LevelOut{1},LevelOut{2},names.states(VarCumSum),names.shocks(ShockstoPlot),namesModels,'',fonts);
    disp(' ')
    disp('------------------------------')
    disp(' ')
end

